/*******************************************************************************

        SSC0143 - Programacao Concorrente
        Prof. Julio Cesar Estrela
       
        T1 - Calculo de PI com 10 000 000 de casas decimais
       
        Grupo:
       
        Daniel de Oliveira Pereira              6427171
        Guilherme Takasawa Yagui                6426698
        Vinicius Marques Stocco                 6427233
       
        Março de 2013
       
********************************************************************************        
                 Algoritmo de Gauss-Lagendre - Versao Concorrente      
*******************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gmp.h>
#include <pthread.h>

// numero de casas decimais corretas de PI
#define PREC_DGT        10000000

// numero de bits minimo para armazenar PI com o numero de casas decimais corretas
#define PREC_BIT        (int) ceil((log(10) / log(2)) * PREC_DGT)

// numero de iteracoes necessario para calcular PI com o numero de casas decimais corretas desejado
#define ITR             (int) ceil(log(PREC_DGT) / log(4))

mpf_t an, an1, an2, bn, bn1, bn2, tn, pn;
mpf_t aux, pi;

/*
 calcula a(n+1)
 a(n+1) = (a(n) + b(n))/2
*/
void *A(void *args)
{     
        mpf_add(an2, an1, bn1);		
	mpf_div_ui(an2, an2, 2);
}

/*
 calcula b(n+1)
 b(n+1) = sqrt(a(n)*b(n))
*/
void *B(void *args)
{     
	mpf_mul(bn2, an1, bn1);
	mpf_sqrt(bn2, bn2);
}

/*
 calcula t(n+1) e p(n+1)
 t(n+1) = t(n) - p(n)*(a(n) - a(n+1))^2
 p(n+1) = 2*p(n)
*/
void *T(void *args)
{
	mpf_sub(aux, an, an1);
	mpf_pow_ui(aux, aux, 2);
	mpf_mul(aux, pn, aux);
	mpf_sub(tn, tn, aux);

	mpf_mul_ui(pn, pn, 2);
}

/*
 calcula PI
*/
void *PI(void *args)
{
	mpf_add(pi, an1, bn1);
	mpf_pow_ui(pi, pi, 2);
	mpf_div(pi, pi, tn);
	mpf_div_ui(pi, pi, 4);
}

void *iterate(void *args)
{
	mpf_set(an, an1);
	mpf_set(bn, bn1);
	mpf_set(an1, an2);
	mpf_set(bn1, bn2);
}


int main(char argc, char *argv[])
{
        clock_t t;              // tempo
	int i;
	pthread_t t_a, t_b, t_t, t_pi, t_iterate; //threads

        // define a precisão desejada para cada variavel
        mpf_set_default_prec(PREC_BIT);
       
        // marca  o clock inicial da execucao do algoritmo
        t = clock();
       
        // aloca memoria para as variaveis utilizadas para calcular PI
        mpf_init(an);
        mpf_init(an1);
        mpf_init(an2);
        mpf_init(bn);
        mpf_init(bn1);
        mpf_init(bn2);
        mpf_init(tn);
        mpf_init(pn);
	mpf_init(aux);
        mpf_init(pi);

        // inicializa a(0) = 1
        mpf_set_d(an2, 1.0);

        // inicializa b(0) = 1/sqrt(2)
	mpf_set_d(bn2, 1.0);
	mpf_set_d(aux, 2.0);
	mpf_sqrt(aux, aux);
	mpf_div(bn2, bn2, aux);

	// inicializa t(0) = 1/4
	mpf_set_d(tn, 1.0);
	mpf_div_ui(tn, tn, 4);

        // inicializa p(0) = 1
	mpf_set_d(pn, 1.0);

	pthread_create(&t_iterate, NULL, &iterate, NULL);
	pthread_create(&t_a, NULL, &A, NULL);
	pthread_create(&t_b, NULL, &B, NULL);
        pthread_join(t_iterate, NULL);
        pthread_join(t_a, NULL);
        pthread_join(t_b, NULL);

        for (i = 0; i < ITR-1; i++)
        {
		pthread_create(&t_iterate, NULL, &iterate, NULL);
        	pthread_join(t_iterate, NULL);
		pthread_create(&t_a, NULL, &A, NULL);
		pthread_create(&t_b, NULL, &B, NULL);
		pthread_create(&t_t, NULL, &T, NULL);
		pthread_create(&t_pi, NULL, &PI, NULL);     
                pthread_join(t_a, NULL);
                pthread_join(t_b, NULL);
                pthread_join(t_t, NULL);           
		pthread_join(t_pi, NULL);
		//gmp_printf("Iteracao %d, pi = %.*Ff \n", i, 100, pi);
        }
       
        // libera a memoria alocada para o calculo de PI
        mpf_clear(an);
        mpf_clear(an1);
        mpf_clear(an2);
        mpf_clear(bn);
        mpf_clear(bn1);
        mpf_clear(bn2);
        mpf_clear(tn);
        mpf_clear(pn);
        mpf_clear(pi);
        mpf_clear(aux);

       
        // calcula o numero clocks necessario para executar o algoritmo
        t = clock() - t;
       
        printf("\nTempo de execucao paralela: %.3f s \n\n", (double) t / CLOCKS_PER_SEC);

        return 0;
}

